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Abstract 

Origin of the magnetic field ubiquitous in the Universe is studied based on the 
Biermann mechanism, which is expected to work in the non-barotropic region. We 
perform a series of two-dimensional M HD simulations of the first generation supernova 
remnant (SNR) expanding in the inhomogeneous interstellar matter (ISM) and study 
the Biermann mechanism working in the interior of the SNR. Especially, we pay 
attention to the relaxation process of electron and ion temperatures via the Coulomb 
interaction. In the early SNR in which the electron temperature is much lower than 
the ion temperature, the Biermann mechanism is ineffective, since the gradient of 
electron pressure is small. Magnetic fields begin to be generated just behind the shock 
front when the electron temperature is sufficiently relaxed. Assuming the explosion 
energy of 10 52 erg, the total magnetic energy generated reaches about 10 26 erg and 
does not depend strongly on the parameters of either SNR or ISM. Analytic expression 
to estimate the magnetic total energy is presented and it is shown this agrees well 
with our numerical results. Finally we evaluate the expected amplitude of magnetic 
fields in protogalaxies as ~ 10~ 19 G, which is sufficient for seed fields of the observed 
galactic magnetic fields. 

Key words: magnetic fields — galaxies: high-redshift — methods: numerical 
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1. Introduction 



The origin of magnetic fields in the universe has been an active field of modern cosmology 
and astrophysics (Widrow, 2002; Kulsrud & Zweibel, 2008). Concerning the origin of galactic 
magnetic fields, which are observationally known to be of order lpG, it is argued that dynamo 
mechanism could amplify and maintain magnetic fields if there were some tiny fields (~ 10~ 20 G) 
at the galaxy formation (for example, Lesch & Chiba, 1995). There have been many approaches 
to the origin of the seed fields such as cosmological scenarios (Turner & Widrow, 1988; Bamba 
& Yokoyama, 2004; Takahashi et al., 2005; Ichiki et al., 2006) and generation via reionization 
(Gnedin et al., 2000; Langer et al., 2005). These scenarios predict that the universe is filled with 
very tiny magnetic fields about 10 -15 — 10~ 25 G and there is a possibility that they could be 
observationally probed through the detection of delayed high-energy emission from gamma-ray 
bursts (Plaga, 1995; Ichiki et al., 2008; Takahashi et al, 2008). 

In this article, we focus on magnetogenesis in primordial supernova remnant (SNR) via 
the Biermann mechanism (Biermann, 1950). The Biermann mechanism was also studied in 
other systems such as the formation of protogalaxies (Kulsrud et al., 1997) and Pop III stars 
(Xu et al., 2008). The Biermann mechanism is an induction process of magnetic fields and can 
be derived by combining a generalised Ohm's law with the Maxwell equation (Widrow, 2002), 

— -Vx(v,xB) = - , (1) 

neglecting the electric resistivity. Here, v p , m e , e, c, P e , and p e are the velocity of protons, the 

mass of electrons, the elementary electric charge, the light velocity, the pressure of electrons, 

and the mass density of electrons, respectively. When we assume n e ~ n p ~ n, P e ~ P/2, and 

v p ~ v where n e , n p , n, P, and v represent the number density of electrons, protons, and ions, 

the total gas pressure, and the fluid velocity, respectively, equation (1) can be rewritten as 

9B „ . „, VP x Vp 

= Vx(vxB) + a 5—^. 2 

at p A 

The coefficient a is defined as a = m p c/2e ~ 0.5 x 10~ 4 G sec, where m p is the mass of protons. 

The Biermann mechanism requires a vorticity of plasma gas. The vorticity is generated 
in the region where the spatial gradients of the pressure and density are not parallel. There 
are some studies of the interaction of the shock of SNR and interstellar clouds (Klein et al., 
1994; Nakamura et al., 2006). These studies revealed that the interaction produces the vorticity 
efficiently in the shocked region. The interaction of inhomogeneous interstellar medium and 
the shock of SNR also drives the turbulent motion of the gas, and the vorticity is generated 
in the bubble of SNR (Balsara et al., 2001). Therefore, primordial SNRs occurred before the 
galaxy formation can be a candidate of the origin of galactic magnetic fields. Furthermore, 
it is important to note that primordial supernova explosions play an important role in the 
magnetogenesis, since the initial mass function (IMF) of population III stars is believed to be 



2 



substantially top-heavy (e.g., Abel et al., 2002). 

For the generation of the seed fields by the primordial supernova explosions, Miranda 
et al. (1998) studied the process with multiple explosion scenario of Population III objects. 
Although their situation considered is more or less similar to ours, their calculations are based 
on several ad hoc assumptions. First, they calculate magnetic field strength assuming a constant 
generation rate which is actually dependent on the configuration of the density and pressure 
gradients. Second, they considered multiple SN explosions, and most of the magnetic fields is 
produced by the hypothetical objects with a mass M = 10 6 - 10 10 M Q , and explosion energy 
E = 10 56 — 10 60 erg. This assumption seems inconsistent with the current understanding of 
star formation. 

In our previous article (Hanayama et al., 2005), we studied the generation of mag- 
netic fields by the Biermann mechanism with realistic two-dimensional magnetohydrodynamic 
(MHD) numerical simulations. A single SNR with the explosion energy 10 53 erg was put in in- 
homogeneous ISM and the generation and evolution of magnetic fields were followed. Then we 
found the total energy of the magnetic fields to be 10 28 ~ 10 31 erg depending on the supernova 
and ISM parameters. 

In this article, we present more detailed analysis especially focusing on the temperature 
relaxation between ions and electrons. When we derived equation (2) from equation (1), we 
used P e ~ P which comes from the assumption that the ions and electrons have the same 
temperature, T e ~ T p . In fact, the temperatures of ions Tj and electrons T e just behind the 
shock differ each other and the Rankine-Hugoniot shock jump condition indicates 

T i = Ar^k, (3) 

16 k-q 
_ 3 m e v^ nh 

Ie ~16 k B ' [V 
where t>bub is the shock velocity. Noting that the time scale of the Coulomb interactions between 
electrons and ions (Spitzer, 1962) is longer than the age of SNR in the early phase, there is a 
possibility that the electron temperature behind the SNR shock front is much lower than the 
ion temperature. In this case, the analysis with equation (2) would overestimate the generation 
of magnetic fields. 

For a long time, the temperature relaxation of ions and electrons behind the shock of 
SNR has been actively discussed by many authors (Itoh, 1978; Cox & Anderson, 1982; Cui & 
Cox, 1992; Masai, 1994). According to Cox & Anderson (1982), in the unequilibrated region 
in the early phase of SNR the electron temperature is spatially isothermal. If this is the case, 
the gradient of electron pressure is negligible and, therefore, the Biermann mechanism is there 
ineffective. Thus, we have to follow the electron temperature in order to analyze the generation 
of magnetic fields correctly. This is the main ingredient of the current article. 

On the other hand, in the radiative cooling phase, although the electron temperature 
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is expected to be sufficiently equilibrated to the ion temperature, the radiative cooling process 
increases its importance with the decrease of the temperature at the shock front, and the 
thermal structure of the SNR becomes isothermal again (Slavin & Cox, 1992). Then, the 
Biermann mechanism does not work effectively in the SNR. Therefore, in this study, we focus 
on the generation of the magnetic fields before the radiative cooling phase. 

In section 2, the basic equations for numerical simulation of primordial SNRs, the nu- 
merical method, and the initial condition of SNR and its surrounding environment are given. In 
section 3, we present the results of numerical simulations. We discuss the analytic estimations 
of the generated magnetic fields in section 4. 

2. Models and Methods 

2.1. The Basic Equations 

The basic equations for our numerical simulations are as follows: 

dp 



dt 



V-(H = 0, (5) 



^ + V.(^) = -VP-v(^) + i(B.V)B, (6) 




+ V 



(E + P)v + — {B x (v x B)} 

4:71 



0, (7) 



VP x Vo 

= V x(,xB) + a^4^, (8) 
at p A 

where p, v, E, P, B, and T are the density, velocity, total energy, pressure, magnetic field, and 

temperature, respectively. Here the total energy of the gas is defined as 

where we fix the value of the adiabatic index 7 to 5/3, which is valid in case of the non- 
relativistic gas. The last term of the right-hand side of equation (8) is called Biermann term 
and expresses the Biermann effect, where a = m p c/2e. We assumed the gas is fully ionized and 
the number ratio of H to He to be 10 : 1. Then the mean molecular weights are p>o — 14/11 = 1.27 
for nuclei and m = 14/23 = 0.61 for fully-ionized ions and electrons. The number density of the 
hydrogen atoms is % = (10/ll)n where n means the number density of the ions (nn + nn e ). 
Since we focus on the generation of the magnetic fields in adiabatic expansion phase, we ignore 
the radiative cooling. 

2.1.1. Generation Region of Magnetic Fields 

We pay attention to the region where the magnetic fields is generated by the Biermann 
mechanism. The interaction between the shock wave and the inhomogeneous ISM generates the 
vorticity. Therefore, especially, in the transition layer of the shock front, the baroclinic term 
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(VP x Vp)/p 2 is expected to work effectively because the gradient of the pressure is much 
larger than in other regions. However, the thickness of the layer is much smaller compared with 
the radius of SNR, and is on the order of, or larger than the ion inertia length of ~ 10 7 /n 1//2 cm 
(Itoh, 1984). In such a region, although the behaviour of the plasma particle must be considered, 
it depends on complex processes of the interaction of ions and electrons, and is conventionally 
studied by particle simulations (Shimada & Hoshino, 2000; Kato & Takabe, 2008). If the 
electron pressure P e is not as high as the proton pressure P p there, the Biermann mechanism 
may not work effectively. In fact, Shimada & Hoshino (2000) predicted the ratio between the 
electron and the ion temperatures T e /Tj ~ 0.2 in the transition layer with Alfvenic Mach number 
of Ma = 20. Observationally, the ratio T e /Ti is less than ~ 0.1 in the transition layer with a 
shock velocity larger than ~ 2000 km s" 1 (Rakowski, 2005). Since in the transition layer the 
electron temperature is expected to be lower than the ion temperature, we ignore the generation 
process working in the transition layer as long as we consider the early phase of SNR. That is, 
we only take the Biermann effect in the interior of SNR into account, even if it would lead to 
an underestimate of magnetic fields. Thus, our analysis is rather conservative. 

In our calculation we introduce a switch C sw that discriminates the region of the 
post shock region where we calculate the Biermann term (a VP x Vp) / p 2 in our one- 
temperature/one-pressure description. For the search of the boundary between the transition 
layer and the post shock region, the switch C sw is defined as follows : C sw = — 1 in the region of 
(VP) ■ v < ; C sw = 1 in the region of (VP) ■ v > ; C sw = in the region of (VP) • v = 0. If 
the switch C sw is equal to or —1, then we reduce the Biermann term to zero. In this way, we 
discuss the generation process of the magnetic fields in the inner region of SNR. Furthermore, 
since a region with T e /T < 0.99 has isothermal T e distribution, magnetic generation is ignor- 
able there. Therefore, we restrict the generation region of the magnetic fields to the postshock 
region with T e /T = g > 0.99. Thus, our calculation must give a same result obtained with 
two-temperature simulation effectively. Then, the induction equation (8) with the relaxation 
effect of the electron temperature is written as 

OB V(oP) x Vp 

— = V x (v x B) + a Ky ' -, (g > 0.99 and C sw = 1) 

8% P " < 10 > 

-j— = V x (v x B). (others) 

2.1.2. Equilibrium Equation of Electron Temperature 

We also consider the relaxation of the electron temperature. In the adiabatic expansion 
phase of SNR, the relaxation equation of electron temperature is derived from the energy 
equation for the electron gas with the heating effect driven by the Coulomb collision with ions. 
Here, we ignore the thermal conduction and the radiative cooling. The relaxation equation was 
given by Itoh (1978), and we convert the formulation according to Cox & Anderson (1982), 

^ = ^ + (-V)/ = ^^j, (11) 



where In A is the Coulomb logarithm and / is a function of g = T e /T defined as 




(12) 



After solving equation (11) for f,g = T e /T can be obtained from an approximate formula 



We solve equation (11) with other MHD equations and obtain the electron temperature T e 
with using equation (13) and the gas temperature distribution T. Since the temperature ratio 
g = T e /T takes 1/1836 just behind the transition region (see equations (3) and (4)), we take a 
boundary condition as g = 1/1836 at the transition region (C sw = —1). 

2.2. Numerical Methods 

We solve the above equations by a two-dimensional MHD code in the cylindrical coor- 
dinates (r,z,(j>) assuming axial symmetry around the symmetry axis (z). Since we start the 
calculation from B = 0, it is sufficient to consider only the 0-component of the magnetic fields 
which is generated from the Biermann term. 

The code is based on the MHD Roe scheme (Roe, 1981; Cargo & Gallice, 1997) coupled 
with MUSCL technique (Hirsch, 1990) to achieve the second order spatial and temporal ac- 
curacy. We employ improved rotated-hybrid Riemann solvers (Nishikawa & Kitamura, 2008) 
to overcome so-called Carbuncle phenomenon, which is a numerical instability appeared in 
the shock wave. Curing V • B error, we adopted the hyperbolic divergence cleaning given by 
Dedner et al. (2002) and Matsumoto (2007). We note that equation (11) is also solved by using 
the fully upwind scheme with MUSCL and Carbuncle cure technique which achieve the second 
order spatial-temporal accuracy and numerical stability. 

The numerical scheme was tested by several problems. We checked the code by com- 
paring with MHD shock tube problems (Brio and Wu, 1988) and the adiabatic SNR with the 
analytic Sedov solution (Sedov, 1959). In high resolution calculation (4000 x 4000 grids) of 
primordial SNR expanding into a uniform ISM (explosion energy E and ISM density n are 
10 51 erg and 0.2 cm -3 , respectively), the postshock density peak value of the shock front at 
1500 yr after the SN explosion was 95.8% of the analytical value (four times larger than the 
density of ISM), for example. In addition, we compared our results with those of Nakamura 
et al. (2006) for the case of interaction of a plane-parallel shock with a cloud to check the 
robustness of our calculation of the vorticity derived from the induction equation (8). As for 
the normalized peak value of total circulation T in the evolution of the cloud distracted by the 
planar shock, it is analytically estimated as —1.77 (Klein et al., 1994) for the model of AS8 
(see Nakamura et al., 2006). The result of Nakamura et al. (2006) is —2.01 while that of our 
calculation is —1.61, and both results are in good agreement with Klein's estimation within 
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In computation, our numerical domain covers a region of 128 pc x 128 pc with 2048 x 
2048 grid points. The grid spacings are Ar = Az = 0.0625 pc and the length scale of density 
fluctuation (8 pc, see below) is resolved with 128 grids. This is sufficient for the convergence 
of the total magnetic energy as we show later. In the simulated region of 128 pc x 128 pc, the 
average density of ISM is approximately constant even in the dark halo structure (Kitayama & 
Yoshida, 2005). 

2.3. Initial Conditions 

We assume inhomogeneous ISM with average number density of no = 0.2 cm -3 for initial 
conditions of our fiducial model, and also consider a model with uq = 0.8 cm -3 as a variation. 
These assumptions are consistent with the previous studies (Kitayama & Yoshida, 2005; Greif 
et al., 2007), in which they obtained initial conditions for Population III SN sites that were 
calculated by self-consistent radiative transfer calculations. The scale length of the density 
fluctuation and the amplitude of the inhomogeneity are not well understood. Recently, Wise & 
Abel (2008) studied the formation process of HII region around a very massive first star with 
3D numerical simulations, and showed the anisotropic expansion and inhomogeneous density 
structure in the HII region with the scale of ~ 1 kpc. In their study, the HII region has relatively 
large density inhomogeneities and the spatial scale of the large inhomogeneity seems to be of the 
order 10 ~ 100 pc. Therefore, we assume here inhomogeneity with relative density fluctuation 
5n = (n — n )/no derived from a Gaussian spectrum with concentration around a scale length 
Ao (32 and 64 pc) and amplitude Aq (0.25 and 0.125), where we assume the power spectrum 
with a peak at the wave number k (= 128 pc/Ao) as A exp[— {k — k Q ) 2 /2ak 2 ]/2irk and set 
the variance = 0.5 to achieve —1 < 5n < 1. Avoiding the biased distribution, a normally- 
distributed random number in the range from 0.9 to 1.1 is multiplied to the amplitude of Sn in 
Fourier space and a uniformly-distributed random number in the range from to 27r is added 
to the phase. Then, the number density is given by using 5n in real space as n = rio(l + Sn). 

For the explosion energy of the supernova, we adopt E = 10 52 erg for the fiducial model 
(Kitayama & Yoshida, 2005; Greif et al., 2007). This explosion energy corresponds to stars 
with mass 200 M which explodes as a pair-instability supernova (Fryer et al., 2001). Also we 
consider a model with E = 2.5 x 10 51 erg as a variation. 

We begin the simulation by adding thermal energy of Eq = 10 52 or 2.5 x 10 51 erg to the 
cells near the origin (r, z) = (0,0) with Gaussian distribution. Initial distribution of thermal 
pressure P SN (r, z) with explosion energy E is defined as 



where <tsn 2 represents variance of Gaussian distribution and we assumed ctsn = 2 pc. 

We set the gas density uniform from the origin (0,0) to 10 pc (= 5o"sn) in the radius, 
assuming the effect of the intense radiation field and stellar wind from a progenitor star. In 
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Table 1. Model parameters. 



Model 


A 


A 


n 


E 




(PC) 




(cm J ) 


(erg) 


A... 


32 


0.25 


0.2 


10 52 


B... 


G4 


0.25 


0.2 


10 52 


C... 


32 


0.125 


0.2 


10 52 


D... 


32 


0.25 


0.8 


10 52 


E. .. 


32 


0.25 


0.2 


2.5 x 10 51 



addition, to connect the homogeneous region and inhomogeneous region smoothly, 5n was 
replaced by 8n' where 

2 

Sn'(x, y,z)=l r + Z 5<7SN ) 6n(r,z), (15) 

in the region from 10 pc (= 5<tsn) to 20 pc (= 10<7sn) in the radius. 
We summarise taken parameters in table 1. 

3. Results 

3.1. Dynamical Evolution of SNR 

First of all, let us give some basic formulae about dynamical evolution of SNR. The 
evolution of SNR is divided into the following three phases, (1) free expansion phase, (2) 
Sedov phase, and (3) radiative cooling phase. The free expansion phase continues until the 
ejecta sweeps up about the same amount of the mass of ejecta M ej in ISM around SNR. The 
transition radius from free expantion to Sedov phase -Rsedov and the transition time isedov are 
written as 

R ^~ 19 H^y l, \mk) m ' (16) 

In the Sedov phase, the expansion of shock wave is well approximated by a self-similar solution, 

i4ub~(^y /5 t 2/5 , (is) 

where -Rbub is the shock radius at time t after the explosion and p is the average mass density 
of ISM. The shock speed is given as the time derivative of -Rbub, 

v h ^ = ^ = 2 -(-) 1/5 t- 3/5 - (19) 



dt 5 \ po 



S 



When the gas at the shock front becomes radiative, a dense shell forms at the outer boundary 
of SNR. The transition into the radiative cooling phase occurs at the time t coo \ and the radius 
.Rcooi where 



n 



"3/ 4 / En \ V8 



'—5.1X10 V^^j [j^j , (20) 

^- 12 °^<ro)" 1/5 (l^) 1/5 ' < 21 > 

where we assumed the radiative cooling is dominated by free-free radiation of H and He (Shull, 
1980). On the other hand, the time scale of the Compton cooling for the primordial SNR is 
given as 

tcom P ~7xl0 6 yr(^) , (22) 

which is larger than t coo i for z ^ 30 in which we are interested. As stated above, in this article, 
we focus on the generation of the magnetic fields in the Sedov phase (tscdov < t < t coo \). 

In fact, there is another important time scale at which the electrons and ions can be 
regarded as a one-temperature fluid. The typical time scale is given by the age at which the 
electron and ion equipartition time just behind the shock front r eq is sufficiently shorter than 
the age, that is, r cq < O.lt (Cox, 1972; Itoh, 1978), 

*-- 3 - 3xl ° 4 ^<o^)" 1/7 (lo^) 3/ " < 23 > 

The mean radius of SNR at this time is written as 

At this time, the mean temperature of electrons is (T e ) = 0.92 (T) where the each mean tem- 
perature is defined with the ion number density n as 

/mWm ^ h n 2 \T e /T s ]R 2 dR . , 

T e T s = ^ R , 25 

V /7 / b n 2 R 2 dR V ; 

_ ^n 2 [T/T s ]R 2 dR 

{)/s ~ J R ^n 2 R 2 dR ' { ' 

where T s represents the gas temperature at the shock front. As already explained, the relaxation 

of the electron temperature is crucial for the Biermann mechanism. 

The above characteristic time scales and radii for the five models adopted here are 

summarised in tables 2 and 3. Here, t e nd in table 2 is the time when we stopped the calculation, 

and -Rend in table 3 is the corresponding radius of SNR. 
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Table 2. Time scale for each epoch. 



1\ InflH 


^bcdov 
Clfl 3 vrl 


L rclax 

HO 4 vr) 


t j 

L cnd 


L cool 
HO 5 vrl 


A 

A. . . 


Q 1 
0. ( 


6.6 


9 1 

0.4 


O.U 


B... 


8.7 


3.3 


3.4 


5.0 


C... 


8.7 


3.3 


3.4 


5.0 


D. .. 


5.5 


1.5 


1.7 


1.8 


E. .. 


17.3 


2.4 


4.1 


4.2 



Table 3. Mean radius for each epoch. 



Model 


-^Scdov 




^cnd 


^cool 




(pc) 


(pc) 


(pc) 


(pc) 


A... 


19 


45 


115 


133 


B... 


19 


45 


115 


133 


C... 


19 


45 


115 


133 


D... 


12 


25 


66 


G7 


E. .. 


19 


30 


93 


94 



3.2. Generated Magnetic Fields 

The evolution of the structure of the gas density and magnetic fields for model A, our 
fiducial model, is shown in figure 1. 

When the blast wave expands into the surrounding inhomogeneous medium and the 
electron temperature is sufficiently relaxed, the magnetogenesis by the Biermann mechanism 
starts to work. For this model, magnetogenesis starts when the SNR expands to the radius 
of ~ 45 pc. At t = 6 x 10 4 yr, the radius of the bubble reaches ~ 60 pc and the anisotropic 
structure is clearly seen, which is induced by the interaction of shock front and the density 
inhomogeneity of the ISM. The amplitude of the magnetic field is ~ 10~ 17 G behind the shock 
front at this time. At t = 2 X 10 5 yr, the shock expands to ~ 90 pc, and the magnetic fields 
are ~ 10 -18 G just behind the shock front while they are 10 -17 G for the inner hot cavity. At 
t = 3.4 x 10 5 yr, the radius reaches ~ 110 pc and the magnetic fields of ~ 10 -17 — 10~ 18 G are 
distributed from the radius of 60 to 110 pc. 

Figure 2 shows the distributions of the gas number density, the electron pressure, the 
magnetic fields and the ratio of the electron temperature T e to the gas temperature T at 
t = 6 x 10 4 yr. At first glance, the distributions of the density (upper- left panel) and the elec- 
tron pressure (upper-right panel) look very similar. However, in fact, the gradient vectors of 
the density and the electron pressure are not exactly parallel (lower panels), which is necessary 
for the Biermann mechanism to work. These panels show that the density gradient has the 
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t =6.06e+04 yr 




t =6.06e+04 yr 



- VP. 



10 20 30 40 50 60 
r[pc] 



10 20 30 40 50 60 
r[pc] 



P 9 [dyn/cm 2 ] 
3.0-1 10 

2.5-1 10 

2.0-1 0'° 

1.5-10" 10 

1.0-10 10 

5.0-10" 



t =6.06e+04 yr 



Vp — VP e -B„[G] 
2-1 0"' 



1-10"' 



t =6.06e+04 yr 




1 -1 1 



2-10 1 




10 20 30 40 50 60 
r[pc] 



10 20 30 40 50 60 
r[pc] 



Fig. 2. Distributions of the gas number density (upper-left), the electron pressure (upper-right), the 
magnetic hclds (lower- left), and the ratio of the electron temperature and the gas temperature (lower-right) 
at t = 6 x 10 4 yr for model A. Red and blue arrows represent directions of the gradients of the density and 
the electron pressure, respectively 

relatively larger azimuthal component compared with the pressure gradient. This would origi- 
nate from the density inhomogeneity in ISM. This is consistent with the picture of Nakamura 
et al. (2006) where it was suggested that the density structure is affected by the distribution of 
the diffuse cloud, while the pressure is mainly determined by the global structure of the SNR 
and has more or less radial gradient. Here it should be important to note that magnetogenesis 
does not occur in the deep interior of the bubble. This is because, as we see in the lower-right 
panel, the electron temperature has not been relaxed due to low density. 

The probability distribution function (PDF) of the magnetic field strength is defined as 

#(log|B,|) 



n(io g |^|) = 



N, 



cells 



(27) 
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Fig. 3. PDF of the strength of the magnetic fields at t = 6 x 10 4 (dashed line), 2 x 10 5 (dotted line), and 
3.4 x 10 5 (solid line) yr for model A. 

where iV(log|S^|) is the number of cells with magnetic fields and iV cc ii s = 2048 2 is the total 
number of cells. The PDFs at these epochs are shown in figure 3. We can see that the magnetic 
fields of ~ 10~ 17 — 10™ 18 G are generated at each epoch. The magnetic fields of ~ 10 -16 G are 
generated in restricted number of cells, although the population is not so large. 

Figure 4 shows the square root of the energy spectrum of the magnetic fields for several 
times. The energy spectrum of the magnetic fields is calculated by using a definition of shell- 
averaged magnetic power spectrum Pm(&) derived by Christensson et al. (2001). Fourier series 
Bf(k r ,k z ) is written as 

B f (k r ,k z ) = S L _J L _ L ^B^z)e-^^\ (28) 

where k r , and k z represent wave numbers for respective directions in a period 2L (= 256 pc). 
We assume B = (0,5^,0) in the first quadrant, and the antisymmetric one against the axis in 
other quadrants. The power spectrum Pm(&) i n Fourier space is defined as 

P M (k) = (B}.B f ), (29) 

where (B*j ■ Bj) is the value averaged over the shells with constant k = \k\. Then, the shell- 
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1 



10 

Mpc] 



100 



Fig. 4. Square root of the energy spectrum of the magnetic fields at t = 6 x 10 4 (dashed line), 2 x 10 5 
(dotted line), and 3.4 x 10 5 (solid line) yr for model A. 

integrated magnetic energy spectrum Euik) is written as 

E M (k) = 2nkP M (k). (30) 

In figure 4, ^Je~m\X) is plotted against the wave length A = 2L/k. At t = 6 x 10 4 yr, the peak 
value is ~ 4 x 10~ 19 G, and the corresponding wave length Abp is 37 pc which is near the scale 
of the initial fluctuation of ISM, Ao = 32 pc. At t = 2 x 10 5 yr and 3.4 x 10 5 yr, the spectra are 
almost identical. The energy spectrum takes the maximum at Abp = 51 pc which is nearly equal 
to 1.5Ao- We can also see that the generation of magnetic fields completed by t = 2 x 10 5 yr. 

The size of the magnetic loop extending in (^-direction equals to the typical coherence 
length appeared in figure 4, ~ Abp/2. Then, the radius of the toroidal fields r t is considered to 
be r t = App/4 ~ 1.5Ao/4, and the coherence length is estimated as 27ir t ~ 1.5Abp ~ 2Ao ~ 64 pc 
at 3.4 x 10 5 yr. This indicates that the magnetic field is expected to be observed with a 
3-dimensional size of 2irr t ~ 64 pc from our 2-dimensional simulations. 

3. 3. Comparison of Models 

In this subsection, we make a comparison of the results of the five models (A-E) shown 
in table 1. The PDFs of the magnetic field strength for the five models are shown in figure 

14 




-19 -18 -17 -16 
log|B,| [G] 

Fig. 5. Comparison of PDFs of the strength of the magnetic fields for each model. Black, blue, green, 
yellow, and red lines represent the models A (t = 3.4 x 10 5 yr), B (t = 3.4 x 10 5 yr), C (t = 3.4 x 10 5 yr), 
D (t = 1.2 x 10 5 yr), and E (t = 4.1 x 10 5 yr), respectively. Dashed (Al) and dotted (A2) lines represent 
t = 6 x 10 4 and 2 x 10 5 yr for model A. 

5. This indicates that the amplitude of magnetic fields are not so different for the variation 
of our models. Considering the maximum amplitude of \B^\ with n > 10~ 2 , 1^(11 > 10" 2 )| 
of model A is roughly twice as large as that of models B and C. The radius and large-scale 
structure for models D (R = 60 pc, t = 1.2 x 10 5 yr) and E (R = 90 pc, t = 4.1 x 10 5 yr) are 
similar to those of model A at t = 6 x 10 4 yr (Al, dashed line; R = 60 pc) and 2 x 10 5 yr (A2, 
dotted line; R = 90 pc), respectively. This simply comes from the self-similar evolution of the 
SNR in the adiabatic phase. Accordingly, we compare the PDF plots of models D and E with 
curves Al and A2, respectively. Then, for model D, the PDF of the magnetic fields can be 
compared with that of model A at t = 6 x 10 4 yr (Al, dashed line). 1-6^(11 > 10 3 ) | of model 
D at t — 1.2 x 10 5 yr is larger than that of model A at t = 6 x 10 4 yr. This difference comes 
from the fact that the post shock pressure of model D at t rc iax = 1.5 x 10 4 yr is larger than that 
of model A at t re iax = 3.3 x 10 4 yr because the amplitude of the generated magnetic fields is 
fundamentally proportional to the pressure gradient. The same argument can be applied to the 
comparison of model E and model A at t = 2 x 10 5 yr (A2, dotted line). For the population of 
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Fig. 6. Comparison of the square roots of the energy spectrum of the magnetic fields for each model. 
Black, blue, green, yellow, and red lines represent models A (t = 3.4 x 10 5 yr), B (t — 3.4 x 10 5 yr), C 
(f = 3.4 x 10 5 yr), D (t = 1.2 x 10 5 yr), and E (t = 4.1 x 10 5 yr), respectively. Dashed (Al) and dotted 
(A2) lines represent t~6x 10 4 and 2 x 10 5 yr for model A. 

II > 10 -2 , the amplitude of model E is less than that of model A at t = 2 x 10 5 yr. This difference 
also comes from the fact that the post shock pressure of model E at t re iax = 2.4 x 10 4 yr is less 
than that of model A at t re iax = 3.3 x 10 4 yr. 

Figure 6 shows the comparison of the square roots of the energy spectrum of the magnetic 
fields for the five models. For models A, C, and E, the curves take their maxima at ~ 51 pc. 
For models B and D, the peak sizes are 85 and 37 pc, respectively. This shows each spectrum 
has a peak near the scale of 1.5Ao except for high-density model D. The peak magnitude of 
model A is 10 -18 G, and that is twice as large as that of models B and C. For high-density 
model D, we can compare the result of model A at t = 6 x 10 4 yr (Al, dashed line), both of 
which have a similar peak size. The peak wavelength of both models is 37 pc, although the 
amplitude of the magnetic fields of model D is twice as large as that of model A (Al). This 
difference seems to come from the fact that the magnetic fields are generated even in the inner 
region for model D since t/t re \ ax is larger compared with other models. This is also seen in the 
comparison of model E at t = 4.1 x 10 5 yr and model A at t = 2 x 10 5 yr (A2, dotted line). The 
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Fig. 7. Time evolution of the total magnetic energy produced by the Biermann mechanism for various 
models. Black, blue, green, yellow, and red lines represent models A, B, C, D, and E, respectively. 

place of the peak wavelength of both models is 51 pc, and the peak energy density of model E 
is slightly lower than that of model A (A2). 

For the coherence length, if the scale of the maximum energy density Abp is equal to 
1.5Ao, the 3-dimensional length of toroidal magnetic field is estimated as 27tAbp/4 ~ 2A . This 
means that the coherence length of the toroidal field of 10 -18 G is estimated as ~ 64 pc for 
models A, C, and E, while it is ~ 128 pc and ~ 32 pc for models B and D, respectively. 

Time evolution of the total magnetic energy is shown in figure 7. The final magnetic 
energy for each model indicated by the figure is, ~ 2 x 10 26 erg for model A, ~ 5 x 10 25 erg 
for model B, ~ 4 x 10 25 erg for model C, ~ 6 x 10 25 erg for model D and ~ 9 x 10 25 erg for 
model E, respectively. Apparent knees around 5 x 10 4 — 2 x 10 5 yr come from the time scale 
of the electron temperature equilibrium shown in table 2. The total magnetic energy of model 
A is several times larger than those of other models. This behavior will be interpreted by an 
analytical estimation of the magnitude of the magnetic fields and the time evolution of the total 
energy in the next section. Although there are many uncertainties in our initial conditions, it 
is implied that the amplitude of the generated magnetic fields does not depend on SNR and 
ISM parameters so strongly. 
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Fig. 8. Convergence study of total magnetic energy E B (left panel) at 3.4 x 10 5 yr. dQ (right panel) 
means the relative convergence error in E B of each model compared with the highest resolution model 
7V gri d = 2048, defined as dQ = \E B (N grid )-E B (N grid = 2Q48)\/ E B (N cl = 2048). 7V grid represents the number 
of grid points in one dimension and the spatial resolution of the simulation. Solid line means the result 
of the models which contain the inhomogeneity (model A), while dashed line means that of the models 
without the density fluctuation (uniform ISM). 

Finally, let us argue the convergence of the results of our simulations. Figure 8 shows 
a convergence study in comparison with the simulations of lower spatial resolution in which 
the numerical box size is taken identically but the grid points are reduced as N gr ^ = 512 and 
1024. Eb and dQ represent the total magnetic energy of SNR at 3.4 x 10 5 yr and the relative 
error in Eb of each low-resolution model compared with that of the highest resolution model 
iV grid = 2048, defined as dQ = \E B (N glid ) -E B {N grid = 2048)]/ E B (N gTid = 2048). In both panels, 
solid line is the result of the calculations of the interaction between the inhomogeneous ISM 
and SNR (model A) and dashed line is that of the evolution of SNR in uniform ISM. Left panel 
could be a measure of the numerical error in magnetic fields of our simulations. We can see 
that, in the highest resolution model, Eb of the model with homogeneous ISM is about 10 5 
times smaller than that of model A and the numerical error of the magnetic fields generated by 
the curvature effect in the interior of SNR is negligible. In right panel, solid line shows that our 
simulation is almost converged and the difference is about 10% for the lower resolution model 
of N grid = 1024. This indicates that roughly 128 grids for the radius of density fluctuation 
r o = Aq/4 = 8 pc are required for a converged calculation. 
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4. Discussion 



4-1. Characteristics of Magnetic Fields 

In this subsection we give an order-of-magnitude estimation of the amplitude of the 
magnetic fields and the time evolution of the total magnetic energy generated by the Biermann 
mechanism. We extend the analysis in Hanayama et al. (2005) considering the time evolution of 
the physical quantities of the SNR bubble and the relaxation of the electron temperature. The 
characteristic pressure at the shock front is given by the ram pressure, P ~ P ram ~ (3/4)po^bub> 
and its gradient is estimated as VP ~ P/L where the pressure scale height can be evaluated 
as a typical shell width L ~ Pbub/10. Noting that the density gradient is determined by the 
fluctuation of ISM, we have Vp ~ A p /(2r ) = 2A p Q /X . Thus, the amplitude of magnetic 
fields generated within a characteristic time, r ~ L/-Ubub, is given by (see equation (8)) 



B 



Bier 



VP x Vp 
a - t 



3v huh A 

~ a 




-1/5 



(31) 

where we used equation (19) and put a ~ 0.5 x 10~ 4 G sec. This is reasonably consistent with 
the value obtained from our numerical simulations. 

On the other hand, the growth of the total magnetic energy during a time interval in 
which an SNR expands from volume V to V + dV can be estimated as 

dE B ~ ^^dV ~ ^i££47rP 2 dP 
8tt 8tt 

. fc^t-i*. (32) 

where dV(= 4irR 2 dR) is a difference in the volume of SNR between two epochs t and t + dt 
and we used equation (18). The total magnetic energy Es{t) contained in a SNR of the age 
t is given by the time integration of dE B (t)/dt from t start to t. Assuming t start ~ t rc i ax f° r the 
sufficient equilibrium of the inner region, we obtain the total magnetic energy as 

* §a 2 A 2 E 

'tstart 

125Ag Po 
9a 2 A 2 E n 

„ (hit- In t s tart 



E B {t) ~ / io - r x dt 
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Fig. 9. The same as figure 7 but added lines of the analytic estimation obtained from equation (33). 
Black, blue, yellow, orange and red lines represent models A, B, C, D, and E, respectively. Solid and 
dotted lines represent the numerical result and the analytic estimation. 

W^=Y ( 33 ) 
Vl0 52 erg/ V 10 / 

See table 2. This analytic estimation is also plotted in figure 9 to compare with our numerical 
simulation. Consistency between these estimations and the numerical results is remarkable, 
explaining not only the qualitative behavior but also the absolute magnitudes. 

4-2. Implication for Seed Magnetic Fields 

Now we estimate the spatially-averaged energy density of the magnetic fields in proto- 
galaxies expected from the first star SNR and consider whether they could be a source of the 
seed fields or not. 

The total number density of the SNe is roughly estimated as 

n SN - (34) 

where p\m, r, and M s are the primordial star formation rate (SFR) of Pop III stars per unit 
volume, the duration of the first star formation and a typical mass of the first stars, respectively. 
As for the SFR of Pop III stars, extrapolating the one by Pello et al. (2004) and Ricotti et al. 
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(2004), we have 



P*,m ~ 6 x 10" 4 M & yr _1 Mpc 




(35) 



where fm is the fraction of the Pop III stars in SFR, and we adopt fm = 0.06 that was 
derived under the assumption that very massive black holes produced from first stars end up in 
supermassive black holes in galactic centers (Schneider et al., 2002). If we assume the formation 
period of the first stars continued from z ~ 20 to 10 (r ~ 0.3 Gyr), the total number density of 
the Pop III SNe can be estimated as 



where the number density is in units of physical scale, not comoving. Primordial star formation 
rate was also estimated by Greif & Bromm (2006). Even taking Pop III and Pop II. 5 of their 
classification into account, we confirmed that our results below do not change so much. 

Taking the typical value of the magnetic energy of a Pop III SNR as that of model 
A, E B ~ 10 26 erg, the spatially-averaged magnetic energy density is estimated as, e# ~ 
10 -42 erg cm -3 . If we assume that galaxies are formed in such a magnetized medium, the 
magnetic energy density in protogalaxies is given by 



where A represents the overdensity of protogalaxies. This means that the average magnitude 
of the magnetic fields becomes B ~ 10 -19 G, which would be enough for the required seed 
field of galactic dynamo (Lesch & Chiba, 1995). Although the coherence length of the order of 
10 — 100 pc estimated here is much smaller than the galactic scale, it would be amplified by 
the galactic dynamo (Poezd et al., 1993; Beck et al., 1994, 1996). It may also be amplified by 
interstellar turbulence dynamo to produce fluctuating components (Balsara et al., 2004). 

Finally let us comment on the three-dimensional effects. In this study, we performed 
two-dimensional MHD simulations. However, because of the assumption of the axisymmetry, 
the generated magnetic fields are restricted to the toroidal component. This makes it rather 
hard to argue the coherence length of magnetic fields. Further, the spectrum of magnetic fields 
would be different in three-dimensional simulations, because the vorticity cascade is different 
in 2D and 3D turbulences. We will present three-dimensional simulations in a separate paper 
but we believe that most of the features of the Biermann mechanism in SNR are captured in 
the present study. 




(36) 



e B A^ 




(37) 
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5. Summary 

In this article, we argued the Biermann mechanism in primordial supernova remnants 
through two-dimensional MHD simulations with the Biermann term. We solved simultaneously 
the relaxation of the electron temperature, which is crucial to the efficiency of the Biermann 
mechanism and was not taken into account in our previous study (Hanayama et al., 2005). 
It was found that magnetic fields begin to be generated from t = t re i ax just behind the shock 
front. The total magnetic energy reaches about 10 26 erg and does not depend strongly on the 
parameters of SNR and ISM. We could understood analytically the dependence of the magnetic 
total energy on the parameters and also the time evolution. Finally we evaluated the expected 
amplitude of magnetic fields in protogalaxies, which would be sufficient for seed fields of the 
observed galactic magnetic fields. 
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